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We derive for Ising spins an off-equilibrium generalization of the fluctuation dissipation theorem, 
which is formally identical to the one previously obtained for soft spins with Langevin dynamics 
[L.F.Cugliandolo, J.Kurchan and G.Parisi, J.Phys.I France 4, 1641 (1994)]. The result is quite gen- 
eral and holds both for dynamics with conserved and non conserved order parameter. On the basis 
of this fluctuation dissipation relation, we construct an efficient numerical algorithm for the com- 
putation of the linear response function without imposing the perturbing field, which is alternative 
to those of Chatelain [J.Phys. A 36, 10739 (2003)] and Ricci-Tersenghi [Phys.Rev.E 68, 065104(R) 
(2003)]. As applications of the new algorithm, we present very accurate data for the linear response 
function of the Ising chain, with conserved and non conserved order parameter dynamics, finding 
that in both cases the structure is the same with a very simple physical interpretation. We also 
compute the integrated response function of the two dimensional Ising model, confirming that it 
obeys scaling x(*>£™) — tw a f{t/tw), with a = 0.26 ± 0.01, as previously found with a different 
method. 

PACS: 05.70.Ln, 75.40.Gb, 05.40.-a 



I. INTRODUCTION 



In the recent big effort devoted to the understanding of systems out of equilibrium, of particular interest is the 
problem of the generalization of the fluctuation dissipation theorem (FDT). The autocorrelation function C(t — t') of 
some local observable and the corresponding linear response function R(t — t') in equilibrium are related by the FDT 

l dC(t-t') 

The question is whether an analogous relation exists also away from equilibrium, namely whether it is still possible to 
connect the response function to properties of the unperturbed dynamics, possibly in the form of correlation functions. 

A positive answer to this question exists when the time evolution is of the Langevin type. Consider a system with 
an order parameter field 4>(x) evolving with the equation of motion 

^M=B(cl>(x,t))+r ] (x,t) (2) 

where B{(f>(x, £)) is the deterministic force and r](x,t) is a white, zero-mean gaussian noise. Then, the linear response 
function is simply given by the correlation function of the order parameter with the noise 

2TR(t,t') = (<j>(x,t)r,(x,t')) (3) 

where T is the temperature of the thermal bath and t > t' by causality. It is straightforward [1] to recast the above 
relation in the form 

where 

A(t,t') = 1 [(<f>(x,t)B(q>(x,t'))) - (B (<f)(x, t)) (p(x, t'))] (5) 

is the so called asymmetry. Eq. (4), or (3), qualifies as an extension of the FDT out of equilibrium, since in the right 
hand side there appear unperturbed correlation functions and, when time translation and time inversion invariance 
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holds, reduces to the equilibrium FDT (1). In Appendix I we show that this equation holds in the same form both 
for conserved order parameter (COP) and non conserved order parameter (NCOP) dynamics. 

The next interesting question is whether one can do the same also in the case of discrete spin variables, where there 
is no stochastic differential equation and, therefore, Eq. (3) is not available. For spin variables governed by a master 
equation, this problem has been considered in recent papers by Chatelain [2], Ricci-Tersenghi [3], Diezemann [4] and 
Crisanti and Ritort [5]. However, although important for computational and analytical calculations, their results, as 
we shall explain below, cannot be regarded as generalizations of the FDT in the sense of Eq. (4). In rcfs. [2,3], a 
scheme is presented where R(t 7 t') is related to unperturbed correlation functions which, however, are not computed 
in the true dynamics of the system. Rather, as will be clarified in Section IV, these correlation functions are computed 
through an auxiliary dynamical rule, instrumental to the construction of an algorithm for the computation of R(t, t 1 ). 
On the other hand, in refs. [4,5], R(t,t') is related to quantities that can be extracted from the dynamics of the 
unperturbed system, but not all of them can be put in the form of correlation functions. 

Instead, we have succeeded in deriving for Ising spin systems a genuine off equilibrium generalization of the FDT, 
which takes exactly the same form as Eqs. (4) and (5) and which holds, as in the Langevin case, for NCOP (spin flip) 
and COP (spin exchange) dynamics. Furthermore, using this result we have derived an efficient numerical method 
for the computation of the response function, without imposing the perturbing conjugate field, which is alternative 
to those of refs. [2,3]. 

The paper is organized as follows: in Section II we introduce the formalism and derive the off equilibrium FDT. 
In Section III we introduce the dynamics in discrete time in order to develop a numerical algorithm based on the 
fluctuation dissipation relation. Section IV is devoted to a comparative discussion of our method with those of 
Chatelain and Ricci-Tersenghi. In Section V the algorithm is applied to the investigation of the scaling properties of 
the response function in the one dimensional Ising model with NCOP (Section V A) and COP dynamics (Section VB). 
As a further application of the method we study in Sec. VI the integrated response function of the Ising model in 
d = 2. In Section VII we make concluding remarks. 



II. FLUCTUATION DISSIPATION RELATION FOR ISING SPINS 



We consider a system of Ising spins Si — ±1 executing a markovian stochastic process. The problem is to compute 
the linear response Rij (t, t 1 ) on the spin at the site i and at the time t, due to an impulse of external field at an 
earlier time t 1 and at the site j. Let 



hj(t) = h5ij6(t - t')9(t' + At-t) 



(6) 



be the magnetic field on the j-th site acting during the time interval [t', t' + At], where 9 is the Heavyside step function. 
The response function then is given by [2,5] 



Ri,j(t,t')= lim — 



1 d( Si (t)} 



At^o At dhj(t') 



h=0 



where 



d(si(t)) 



dh 3 {t>) 



h=0 



v rr~i fih'i f i Art dp h ([s'},t' + At\[s"},t') 

2^ Sip{[s\,t\[s),t +At) — 



p{\s"W) 



(7) 



(8) 



h=0 



and [s] are spin configurations. 

Let us concentrate on the factor containing the conditional probability in the presence of the external field p h ( [s'] , t' + 
At\[s"], t'). In general, the conditional probability for At sufficiently small is given by 

p([s],t + At\[s'],t) = 5 [sUsl] + w([s'} -> [s])At + <D(At 2 ), (9) 

where we have used the boundary condition p([s], t\ [s'], t) = 8[ s \,[s'\- Normalization of the probability implies 

5>(M -[*']) = <). (io) 



Furthermore, the transition rates must verify detailed balance 

w([s] [s'])exp(-W[ S ]/T) = w([s'} [s])exp(-H[s'}/T), 



(11) 
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where H[s] is the Hamiltonian of the system. In the following we separate explicitly the diagonal from the off-diagonal 
contributions in w([s] — ► [s']) 

*>([*] ^ [*']) = -*[,],[.>] £ #]-K']) + (i-^],mMW-M), (12) 
[,»]*{,] 

where we have used Eq. (10). 

Introducing the perturbing field as an extra term AH[s] = —Sjhj in the Hamiltonian, to linear order in h the most 
general form of the perturbed transition rates w h ([s] — > [s']) compatible with the detailed balance condition is (see 
Appendix II) 

w h ([s] [a']) = w°([s] Is'}) {l - ^(s, - s'j) + M([a], [,'])} , (13) 

where M([s],[s']) is an arbitrary function of order h/T symmetric with respect to the exchange [s] <-> [s'], and 
— * [s']) are unspecified unperturbed transition rates, which satisfy detailed balance. Notice that, since Eq. (11) 
reduces to an identity for [s] ^ [s'], Eq. (13) does not hold for the diagonal contribution w' t ([s] — > [s]) which, in turn, 
can be obtained by the normalization condition J2[ s '] wh ([ s ] ~ * I s ']) = 0- In the following, for simplicity, we will take 
M([s], [s 1 ]) = and the role of a different choice for M([s], [s']) will be discussed in sec. IV. 
Using Eqs. (9), (12) and (13) we obtain 



T dp^([s],t + At\[s'],t) 



dhj 



h=0 [«"]#[«] 



(14) 

and inserting this result in Eq. (8), the response function can be written as the sum of two contributions [6,7] 

TR itj (t, t') = ^im [TD it j(t, t', At) + TDi t j(t, t', At)] , (15) 

where the first term comes from the diagonal part of Eq. (14) 

TD iJ (t,t',At) = \ s iP ([s},t\[s'},t' + At) E W°(M^M)«-<)P([«V) (16) 
W.M [«"M«'] 

whereas Djj takes all the off-diagonal contributions 

TD^{t,t',At) = \ J2 s iP ([s},t\[s'},t' + At)(^ - - [*>([*"], *0- (17) 

W.M.[«"]?*[«'] 

Using the time translational invariance of the conditional probability p([s], t\ [s'], t' + At) = p([s],t — At\[s'],t'), one 
can write Dij(t, t' , At) in the form of a correlation function 

TD itj (t,t',At) = —(siit-A^Bjit')), (18) 

where 

Bi = -£te-*>>M [*"])■ (19) 

[«"] 

Using Eqs. (9) and (12) the off-diagonal contribution can be written as 

TD i At,t',At)= 1 - ACi ^ t ' ) (20) 

where 
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AC<,,(M') - (si(t)[ Sj (t> + At) - Sj (t')}) = 

- s'j)p{[s],t\[s'],t' + At)p([s'},t' + At\[s"],t')p([s"},t'). (21) 



s , s' , s" 



Therefore, putting together Eqs. (18), (20) and taking the limit At — > we obtain 

TRijW) = \ 9Ci '$ ,n \{ S i{t)B 3 {t>)). (22) 
In order to bring this into the form of Eqs. (4) and (5), we notice that from Eqs. (9, 12) follows 

d(sj(t)) _ v dp([s],t) 
dt ^ 3 dt 

[«] 

= -£ £ £ «X([»1-W)p(M,t). (23) 

W [«"]*[«] [«].[«']*[«] 
Hence, after the change of variables [s] — > [s"], [s'\ — > [s] in the second sum, one obtains 

= -£ £ (sj *>>] - K]MW,t) = (BjW). (24) 

M [»"]*[«] 

In a similar way, it is straightforward to derive 

at 

and subtracting this from Eq. (22) we finally find 

TRij(t,t)-- ^ Wt A itj (t,t) 

where Aij (t,t') is given by 



(B i (t)s j (t'))=0 (25) 



(26) 



A Uj (t,t') = - [{sMBtf)) - (B^stf))] . (27) 



Eqs. (26) and (27) are the main result of this paper. They are identical to Eqs. (4) and (5) for Langevin dynamics, 
since the observable B entering in the asymmetries (5) and (27) plays the same role in the two cases. In fact, Eq. (24) 
is the analogous of 

= <*(*(*,*))> (28) 

obtained from Eq. (2) after averaging over the noise. 

In summary, Eq. (26) is a relation between the response function and correlation functions of the unperturbed 
kinetics, which generalizes the FDT. Furthermore, Eq. (26) applies to a wide class of systems. Besides being obeyed 
by soft and hard spins, it holds both for COP and NCOP dynamics. Moreover, as it is clear by its derivation, Eq. (26) 
does not require any particular assumption on the hamiltonian nor on the form of the unperturbed transition rates. 



III. DYNAMICS IN DISCRETE TIME: THE NUMERICAL ALGORITHM 



We now discuss the numerical implementation of the fluctuation dissipation relation derived above, as a method 
to compute the response function without imposing the external magnetic field (6). Let us recall that Eq. (22) was 
obtained letting At — > 



TR, j(t,t') = I lim 



AC,, 3 (M') 
At 



(siit-A^Bjit')) 



(29) 
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In the simulations of an TV-spin system, time is discretized by the elementary spin updates. Measuring time in 
montecarlo steps, the smallest available time e — l/N is the one associated to a single update. Then, in discrete time 
Eq. (29) reads 

TRijM = I QiM+lLzQ&n l {Sl(t e) Bj (t>)) (30) 

and we use this for the numerical calculation of the response function. For completeness we also give the expression 
for the integrated response function 

Xi,j{t,[tM) = J t RiAt't'W (31) 

which correspond to the application of a constant field between the times t w and t. From Eq. (30) we have 

T X ij{t, [t,t«J) = TeY, RijW) = 2 [<7 « ( *'* ) - C ^t w )] - \ E ^ - e )^(*')> ( 32 ) 
t'=t w t'=t w 

where t — e >t > t w > 0, and Y?t'=t m stan ds for the sum over the discrete times in the interval [in,,*]. 



IV. COMPARISON WITH DIFFERENT ALGORITHMS 



Expressions for the response function in a discretized time dynamics have been derived previously by Chatelain [2] 
and Ricci-Tersenghi [3] . Restricting to transition rates of the heat-bath form and to the case of single flip dynamics 
they have obtained 

t 

Txi,j(t, \t,t w ]) = Yl Ewi^wfew-^fw))^). (33) 

I(t,t w ) r=t ™ 

where I(t) is the index of the site updated at the discrete time r, I(t, t w ) is a specific sequence of /(r)'s between the 
times t w and t, sj v (r) = tanh[/3/i]^(r)] and (t) is the local field due to the spins interacting with Sj. 

It is important to stress the differences between Eqs. (33) and (32). Although in the r.h.s. of Eq. (33) there 
appears an unperturbed correlation function, this is computed in the ad hoc kinetic rule introduced for the purpose of 
evaluating the response function. In Eq. (32), instead, the average in the r.h.s. is computed in the true unperturbed 
dynamics of the system. This important difference arises because of the presence of the delta function 5j,j( r ) in 
Eq. (33), which constrains to update at the time r only the spin at the site j, where the external field is applied. 
Then, in the averaging procedure, only the subset of dynamical trajectories with I(t) — j are considered, while all the 
others get zero statistical weight, which is not what happens in the true unperturbed dynamics. Eq. (33), therefore, 
although useful for the computation of the response function, is operatively restricted to a numerical protocol with 
a sequential updating satisfying the constraint imposed by the delta function Sjj( T )- The correlations functions 
appearing in this equation cannot be extracted numerically from the behavior of the original unperturbed system, 
and cannot be accessed in an experiment. 

Another difference between Eq. (33) and Eq. (32) concerns the choice of M([s], [s']). Our results are obtained with 
M([s], [s']) — 0. Instead, Eq. (33) corresponds to M([s], [s']) ^ 0. In fact, Eq. (33) assumes heat bath transition rates 
w (i s ] ~^ I s ']) = { cx p[~ ^[ s ']/Tlj7{ ex P[ — + ex P[~ ^[ s ']/^]}i expanding this expression to first order in powers 
of h/T, and comparing with Eq. (13), one has 

h. H[s']/T _ H[s]/T 

M([s], [s']) = ^(4 - S3 )- HWW — mjT . (34) 
Retaining M([s], [s']) in Eq. (13) and following the same steps as for M([s], [s']) = 0, one finds the extra term 

AR(t,t>) = L s iP ([s},t\[slt')M([s'}, [*"]) [«,>"] - H)P([*'V) - - W'})p(W},t>)] (35) 

3 r.i r.'i r»«i 
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in addition to the quantities already present on the r.h.s. of Eq. (26). This term cannot be related to correlation 
functions. It is generally believed that the large scale - long time properties of the dynamics (perturbed or not) do 
not depend too much, within a given universality class, on the form of the transition rates. Then one expects the 
corrections AR(t,t') introduced by different choices of M([s], [s']) to be negligible. Indeed, as will be shown in the 
following sections V,VI, numerical results obtained for the Ising model in d = 1,2 with the two algorithms arc not 
sensitive to the choice of M([s], \s'\) . 

V. RESPONSE FUNCTION OF THE ISING CHAIN 

As an application of the numerical method, we compute the autoresponse function R(t, t') = Ri,i(t, t') in the d = 1 
Ising model, with and without conservation of the order parameter. We consider the system prepared in the infinite 
temperature equilibrium state and quenched to the finite temperature T > at the time t = 0. Since the critical 
temperature vanishes for d = 1, the final correlation length £ eq is finite and equilibrium is reached in a finite time 
T eq ~ £| , where z is the dynamic exponent. For deep quenches r eq is large and, after a characteristic time t sc , a well 
defined non-equilibrium scaling regime sets in for t sc < t < r eq , characterized by the growth of the domains with a 
typical size L(t) <~ t x / z . We study the scaling properties of the response function R(t, t') when both t and t' belong 
to the scaling regime. 

A. Non conserved dynamics 

The linear response function in the model with single spin flip dynamics has been computed analytically [8-10]. 
This case, therefore, is useful as a test for the accuracy of the algorithm. In the aging regime t' < t < r eq the 
autoresponse function R(t, t') is given by [9] 

TR(t, t') = e-(*-*')/o (t - t') e- 2t ' [7 (2f) + h {2t')\ (36) 

where I n {x) are the modified Bessel functions. 

In order to improve the signal to noise ratio, we have extracted R(t, t') from the integrated autoresponse function 
X(t, [f + S, t']) by choosing 5 in the following way. Expanding for small 6 we have 

X (t,[t' + S,t'}) 5 dR(t,t>) 
5 " R{t > 1 ] + 2^F~ 

then, for a given level of accuracy, R(t,t') can be obtained from x(t, [t' + S,t'])/8 if 6 is chosen appropriately small. 
Notice that, assuming scaling R(t,t') = t~( a+1 ^ f{t' /t), from Eq. (37) one has that, for a given value of x = t'/t, the 
second term on the r.h.s. produces a relative correction AR(t,t')/R(t,t') = (l/2)[f'(x)/f(x)](S/t) of order S/t. In 
our simulations we have chosen 8=1 and, since the simulated times are t > 100, we have always 6/t < 10 -2 . In the 
following it is understood that all numerical results for R(t, t') are obtained in this way. 

In Fig. 1 we compare the numerical results obtained by means of the algorithm (32), for three different values of t', 
with the exact solution (36) . We have also plotted the data obtained using the algorithm (33) of Ricci-Terscnghi [3] , 
finding an excellent agreement between the curves generated by the different algorithms and the analytical expres- 
sion (36). 

The physical meaning of the exact solution can be understood by replacing Eq. (36) with the simple interpolating 
formula 

TR(t, t') = Azt'- 1 ' 7 -^ -t' + to) 1 / 2 - 1 (38) 

obtained by replacing the Bessel functions with the dominant term in the asymptotic expansion and by inserting to 
as a regularization of R(t,t') at equal times. For NCOP z = 2. With A 2 = l/(y/2n) and to = l/(27r) the simple 
algebraic form (38) gives a very good approximation of the exact solution all over the time domain, from short to 
large time separations. Rewriting Eq. (38) as 

TR(t, t') ~ L(t'y 1 TR sing (t - t') (39) 

where 

TR sing {t-t') = A z {t-t') 1 ' z - 1 (40) 
6 



(37) 



the physical meaning becomes clear, since L(t') -1 is proportional to the density of defects pit') at time t', and 
Rsing(t — t') can be interpreted as the response associated to a single defect. In other words, Eq. (39) means that 
the total response is given by the contribution of a single defect times the density of defects at the time t' . As a 
matter of fact, Eqs. (39) and (40) are the particular realization of a general pattern for the aging part of the response 
function [11] 

R(t, t') ~ Lit'^Rsingit - t')fit/t'). (41) 

The presence of the scaling function f(t/t') in Eq. (41) for d > 1 can be explained as follows: in d = 1 interfaces are 
point like and the interaction between them always produces annihilation. This is accounted for by the defect density 
Lit') -1 , so f(t/t') = 1. In higher dimension, however, defects are extended objects whose interaction can produce a 
wealth of different situations, which are globally described by a suitable scaling function f(t/t'). 

B. Conserved dynamics 

The generality of the structure of Eqs. (39) and (40) may soon be tested by looking at the response function in 
the Ising chain with COP dynamics. While for NCOP the system enters the scaling regime almost immediately, 
since t sc ~ 1, for COP the time t sc for the onset of the scaling regime is of the order of the characteristic time 
T ev ~ exp(4J/T) for the separation (evaporation) of a spin from the boundary of a domain. The equilibration 
time, instead, is given by T eq ~ exp(10J/T) [12]. In order to have a large scaling regime, namely T eq 3> t ac , it is 
necessary to take T/J <C 1, and to choose t w > t sc . Simulations of the system in these conditions are excessively 
time demanding with a conventional montecarlo algorithm. Therefore we have resorted to the fast algorithm of Bortz, 
Kalos and Lebowitz [13], which is much more efficient at low temperatures. With a conventional algorithm a number 
of attempts proportional to r ev is necessary on average before the evaporation of a spin from a domain occurs. Then, 
at low temperatures, a huge amount of attempted moves are rejected, causing a very low efficiency. The algorithm of 
Bortz, Kalos and Lebowitz, instead, is rejection free: moves are always accepted and time is increased proportionally 
to the inverse probability associated with them. We stress that this is not an approximate kinetics, but a clever 
implementation of the exact dynamics. 

From the unperturbed system the response function is extracted through Eq. (32), as in the case of NCOP. About 
the choice of S, for COP the simulated times are t > 10 7 and we can have S/t < 10~ 2 , as required for the correction 
term in Eq. (37) to be negligible, with S — 10 5 . We have performed simulations with T — 0.3 J, corresponding to 
t sc — 6 • 10 5 and r eq ~ 3 • 10 14 . In this conditions the scaling regime is very large. Actually, after a very narrow initial 
regime, where single spins diffuse until they are adsorbed on an interface, no evaporations occur and the system is 
frozen up to times of order r ev . In this time regime the density of defects p(t) stays constant, as shown in Fig. 2. 
Then, for t > T ev , the evaporation-condensation mechanism takes place and the systems gradually enters the scaling 
regime. The range of times explored for the computation of the response function is shown in Fig. 2. This has been 
chosen as a compromise between the necessity to go to the largest accessible times, in order to work well inside the 
scaling regime, and to speed up the simulation to have a good statistics. For COP, the observation of the asymptotic 
behavior p(t) <~ t 1 / 2 , with z = 3, requires very large time [14]. In the range of times explored for the computation of 
Rit,t') the effective exponent z e ff = — [d\og pit) / dlogt]^ 1 has reached the value z e ff = 3.44. 

We have plotted in Fig. 3 the numerical data for Rit, t') together with the curves obtained from the analytical 
form (38), where we have substituted for z the above value of z e ff and we have used A3 and to as fitting parameters. 
The comparison is good and suggests that the physical interpretation, behind the form (39) and (40) of the response 
function, applies also to the d = 1 Ising model with spin exchange dynamics [15]. We expect the more general 
form (41) to hold in higher dimension with COP. 

VI. RESPONSE FUNCTION OF THE D = 2 ISING MODEL 

As a further application of the numerical method, we compute the zero field cooled magnetization (ZFC) %(t, t w ) = 
X(t, [t, t w ]) in the d = 2 Ising model with NCOP, quenched from the infinite temperature equilibrium state to a 
temperature below T c . This quantity has already been measured both by applying the perturbation [11,16-18] or 
by means of the algorithm of Ricci-Tersenghi [3] . In Fig. 4 we compare results obtained with our method and with 
that of Ricci-Tersenghi, for several values of t w in the scaling regime. The agreement between the two algorithms is 
excellent also in this case. The equivalence of the two algorithms both in d = 1 and in d = 2 suggests, recalling the 
discussion at the end of Sec. IV, that different choices of M([s], [s']) do not produce significant differences. 
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Let us comment on the behavior of x(<,t UI ). As it is well known, in the late stage of phase-ordering the interior of 
the growing domains is equilibrated, while interfaces are out of equilibrium. Then, a distinction can be made between 
bulk and interface fluctuations. Accordingly, for the ZFC one has [21-23] 



X(t,t w ) = Xst{t,t w ) + Xag(t,t w ). 



(42) 



Here Xst{t,t w ) is the contribution from the bulk of domains, which behaves as the equilibrium response \eq{t,t w ) in 
the ordered state at the temperature T. This quantity, starting from zero at t = t w , saturates to the value 1 — M 2 , M 
being the equilibrium magnetization. The other term appearing in Eq. (42) , namely the additional aging contribution 
due to the interfaces, is much less known. It is expected to scale as 



In previous studies [11,17,18], an auxiliary dynamics, which prevents flips in the bulk, was used in order to extract 
the aging part of the response in Eq. (42). Here, instead, we have chosen a different technique to isolate Xag(t,t w ): 
We compute the full x(i, t w ) in the Glauber dynamics working at a sufficiently low temperature where \st{t,t w ) is 
negligible. In fact, by choosing T = J, the asymptotic value of Xst(t, t w ) is 1 — M 2 ~ 0.0014, much smaller than the 
computed values of x(t, t w ) in the range of times considered. Then one has t w ) ~ Xag(t, t w ). 

Besides this difference, previous results [11,17,18] on Xag(t,t w ) were obtained with the usual method where a 
perturbation is applied. The strength h of the perturbation must be chosen sufficiently small to work in the linear 
regime. However, by reducing h the signal to noise ratio lowers, and the results get worst. Then one usually runs a 
series of preliminary simulations in order to determine the largest value of h compatible with the requisite of working 
in the linear regime. While this point may be subtle, in the result presented in this paper the limit h — > is taken 
analytically in the derivation of the algorithm. 

We have extracted a from the data of Fig. 4 by plotting x(i, t w ) against t w with x — t/t w held fixed. The results 
are shown on a double logarithmic plot in Fig. 5. According to the scaling form (43), for different values of x the data 
must align on straight lines with the same slope a. This is very well compatible with the curves of Fig. 5, indicating 
that scaling is obeyed. Computing a as the slope of these curves we find a = 0.26 ± 0.01. This result agree very 
well with the value found in [11,17,18]. Once this exponent is known, one obtains data collapse by plotting t^x(t,t„,) 
against x, as shown in the inset of Fig. 5, confirming the validity of the scaling form (43). 



In this paper we have derived a generalization of the FDT out of equilibrium for systems of Ising spins, which 
takes exactly the same form of the FDT generalization previously derived [1] for soft spins, evolving with Langevin 
dynamics. 

We have shown that this fluctuation dissipation relation, which reduces to the usual FDT when equilibrium is 
reached, is obeyed in complete generality both by systems with COP and NCOP. In addition to the theoretical 
interest, as a contribution to the understanding of the FDT in the out of equilibrium regime, our result is promising 
also as a convenient tool for the computation of the linear response function in numerical simulations without applying 
the perturbation, along the line of refs. [2,3]. With standard methods, the requirement to work in the linear regime, 
namely with an adequately small perturbation, sometimes is very subtle and hard to be checked. This problem is 
avoided by this new class of algorithms. Moreover, the statistical accuracy of the results is, for comparable cpu times, 
much better because simulations of perturbed systems usually require additional statistical averages over realizations 
of the (random) perturbation. We have demonstrated the high quality of the results produced by our algorithm by 
computing the response function of the Ising model in d = 1 and the integrated response function in d = 2. In 
d = 2 our results agree with those obtained with the algorithm of ref. [3] and with previous simulations performed 
applying the perturbation. We confirm that x{t->tw) obeys a scaling form (43) with a = 0.26 ± 0.01, in agreement 
with previous determinations [11,17,18] of this exponent. In d = 1, for NCOP our results are in excellent agreement 
with the exact analytical solution and with the simulation made with the algorithm of ref. [3]. In the case of COP, 
where no analytical solution is available, we have obtained results which substantiate the existence of the common 
structure (39,40) of the response function for COP and NCOP. These results show that the algorithm is efficient 
enough to give access to the direct measurement of the impulsive response function R(t,t'), which is too noisy to 
be computed with standard methods. For this reason, previous numerical studies [11,16,19,20] have been necessarily 
directed to the investigation of the integrated response functions, such as the thermoremanent magnetization or the 
zero field cooled magnetization, which are easier to compute. However, as discussed in detail in [11], it is quite delicate 



Xag(^:^«i) — t w /( ). 



(43) 
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a task to extract the properties of R(t, t') from those of the integrated response functions. Therefore, the feasibility 
of a direct computations of R(t, t') is an important development in the field, which is expected to solve a number of 
problems still open [11] on the scaling behavior of R(t, t') for d > 1. 



ACKNOWLEDGMENTS 

We are much indebted to Claudio Castellano for valuable suggestions on the numerical techniques. This work has 
been partially supported from MURST through PRIN-2002. 



APPENDIX I 



Let us write the Langevin equation in the general form 

d(j){x, t) 



dt 



= (zV) p [B(<f>{S,t))+h{sH,t)\ +v{x,t) 



(44) 



where h(x, t) is the external field conjugated to the order parameter, p = or p = 2 for NCOP or COP, respectively, 
and the noise correlator is given by 



(r?(f, t) v (x', t')) = (iV) p 2T<5(£- x')5(t - t'). 
Fourier transforming with respect to space, these become 

d<p(k, t) 



dt 



= k p 



B ([4>}, k,tj +/»(*.*)] +V(k,t) 



and 



{r](k, t)r](k', t')) = 2Tk p (2Tr) d S(k + k')S(t - t') 



(45) 



(46) 



(47) 



where B ([4>], k, tj is the Fourier transform of B (4>(x, t)). 
The linear response function is defined by 



R(k,t,k',t' 



5{4>{k,t)) h 



Sh(k',t') 



(48) 



h=0 



with t > t' . Notice that, since rj{k,t) and k p h(k,t) enter the equation of motion (46) in the same way, we have 

/ S<f>(k,t) \ _ 1 6(<Kk,t)) h 



\5ri{k',t')l k'P 5h{k',t') 



h=0 



where (•) denotes averages in absence of the external field. Then, using the identity [24] 

{ ^k,t)r 1 {k',t')) = 2Tk' p l^^\ 

we find 

2TR(k,t,k',t') = (<f)(k,t)r](k',t')} 

or in real space 

2TR(x,t,x',t') = (<t>(x,t)ri(x' ,t')) 
showing that Eq. (3) holds in the same form for NCOP and COP. 



(49) 

(50) 

(51) 
(52) 
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APPENDIX II 



The transition rates, with and without the external field, satisfy the detailed balance condition (11). Writing 
w h ([s] -> [s']) = w°([s] -> [s'}) Aw{[s] -» Is'}) and H h [s] = H°[s] + AH[s], from Eq. (11) follows 



Aw([s] -> [s'])exp 



AH[s] 



Aw([s'] -> [s])exp 



AW [a'] 



(53) 



Using AH[s] = -Sj/ij, Eq. (53) is satisfied by Aw([s] -> [s']) = cxp [-(l/(2T))/i j (s j - s$)] up to a factor M([s], [s'}) 
which satisfies A4([s], [s']) = M([s'} 7 [s]). Therefore, the most general form of the perturbed transition rates, compat- 
ible with detailed balance, is given by 



Aw([s] -> [s']) = cxp 



M([s],[s'}). 



(54) 



For h — the condition Aw([s] — > [s']) = 1 implies -M([s], [s'}) = 1. Therefore, to linear order in the perturbation 
one has A4([s], [s']) ~ 1 + M([s], [s']). Inserting this result in Eq. (54), and expanding also the exponential term, to 
linear order in h/T one obtains Eq. (13). 
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10 , 

10 10 10 

t (mcs) 

FIG. 2. p(t) in the d = 1 Ising model with COP, T = 0.3 J and J = 1. The number of spins is N = 10 4 . The range of time 
used in the simulations for the computation of the response function is in between the vertical lines. The dashed line represents 
the asymptotic law p(t) ~ t^ 1 ^ 3 . 
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FIG. 4. \(t,t w ) in the d = 2 Ising model with NCOP, T = J = 1. The number of spins is TV = 1600 2 , 
t w = 1 • 10 3 ,t w = 1.5 • W' i ,t w — 2 ■ 10 3 ,t w = 2.5 • 10 3 ,t w = 3 • 10 3 , from top to bottom. Circles represent data obtained 
with the algorithm of Eq. (32), continuous lines are the results with the Ricci-Tersenghi method (33). 
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